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Abstract 

One dimensional propagation of ultraslow optical pulses in an atomic Bose-Einstein conden- 
sate taking into account the dispersion and the spatial inhomogeneity is investigated. Analytical 
and semi-analytical solutions of the dispersive inhomogeneous wave equation modeling the ultra- 
slow pulse propagation are developed and compared against the standard wave equation solvers 
based upon Cranck-Nicholson and pseudo-spectral methods. The role of curvature of the trapping 
potential of the condensate on the amount of dispersion of the ultraslow pulse is pointed out. 
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I. INTRODUCTION 



Possibility of dramatic reduction of speed of light, as demonstrated in an atomic Bose- 
Einstein condensate (BEC) jl| using electromagnetically induced transparency (EIT) 3l, 
opens the way towards practical realization of quantum optical information storage |3|, |j, 
^ y, [7]. In BECs, ultraslow group velocities allow for storage times of coherent optical 
information of the order of few microseconds. To enhance the storage capacity and duration 
it is necessary to inject multiple pulses to be simultaneously present in the condensate. Major 
limiting factor against this goal is the narrow EIT window and group velocity dispersion of 
the slow pulses. 

Propagation of the ultraslow wave packet in one dimensional inhomogeneous dispersive 



atomic condensate can be described in terms of the slowly varying pulse envelope [8| 

dE , I BE ^ , ^d'^E 

where a{z) is the pulse attenuation factor; Vg{z) is the group velocity, and h2{z) is the 
group velocity dispersion. The third order dispersion is found to be much smaller and 
neglected Previous studies are based upon either analytical solutions obtained for an 
effective homogeneous region at the center of the condensate where the dispersion is highest 
or numerical propagation of the pulse with the standard methods such as Cranck-Nicolson 
or pseudo-spectral schemes {9]. 

Present work aims to solve Eq. ([1]) analytically so that optimum experimental conditions 
can be efficiently determined to reduce the dispersive effects to enhance coherent information 
storage capacity of the condensate. In addition to exact analytical expressions of the pulse 
envelopes, an approximate method of determination of the pulse envelope by polynomial or 
Gaussian function fitting to condensate density profile is introduced. The density profile 
only contributes as an integrand to the optical parameters so that the fitting can be done 
with very high accuracy for a simple and fast evaluation of the required integrals. Besides, 
the spatial dependence of the coefficients in the wave function is due to the inhomogeneous 
condensate profile for which non-condensate, thermal component makes a small contribution 
at low temperatures. The practical semi-analytical approach outlined as above is tested 
against Cranck-Nicolson and pseudo-spectral wave equation solvers and shown to be highly 
efficient. The approach is employed to investigate the effect of the curvature of the trapping 
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potential of the condensate, translated to the curvature of the dielectric function of the 
condensate, on the amount of the dispersion of the optical pulse. 

The paper is organized as follows: In Sec. II, optical properties of an atomic BEC system 
under EIT conditions is briefly reviewed. In Sec. Ill, analytical and semi-analytical solutions 
are developed for the wave equation. Numerical schemes used against the semi-analytical 
solution are described in Sec. IV. Results are discussed in Sec.V. Finally, we conclude in 
Section VI. 



II. OPTICAL PROPERTIES OF AN ATOMIC BEC UNDER EIT CONDITIONS 

Atomic condensates can be considered as two component objects, composed of a conden- 
sate cloud in a thermal background so that one can write its density being 

Pir) = pc{r) + prir), (2) 

where jl^ 

Pcir-) = ^^^^e(/i - V{r))e{T^ - T), (3) 



PT{r) = ^^^3 . (4) 
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Here Uq = Anh as/m, m is atomic mass, and is the atomic s-wave scattering length. B(.) 
is the Heaviside step function, gn{x) = Sj-x-^/j", Xt is the thermal de Broglie wavelength, 
(3 = l/ksT, and Tq is the critical temperature. We assume an external trapping potential 
in the form V{f) = {l/2)m{uj'^r'^ + uj'^z'^) with the radial trap frequency and Uz the 
axial trap frequency in the z direction, fi is the chemical potential. At temperatures below 
Tc, the chemical potential /i is determined by yu(T) = pTF{No/Ny^^, where [xtf is the 
chemical potential evaluated under Thomas-Fermi approximation. The condensate fraction 
is evaluated by [hJ No/N = 1-x^- sC(2)/C(3)a;2(l - x^)^/^ with x = T/T,, and C is the 
Riemann-Zeta function. The scaling parameter s is given by s = fiTp/kBTc- 

EIT susceptibility for an atomic BEC of atomic density p can be expressed as 
X{r,uj) = p{r}xii^) with 

( ) = i(r2/2-iA) 

> eoh [(12/2 - iA)(r3/2 - iA) + ni/A] ' ^ ^ 
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where A = u — uq is the detuning of the probe field frequency uj from the atomic resonance 
ujq. Qc is the Rabi frequency of the control field; 7 is the dipole matrix element for the probe 
transition. r2 and denote the dephasing rates of the atomic coherence. 

Through the atomic density, the optical response of the atomic condensate becomes spa- 



tial 
by 



y inhomogeneous, so does the parameters in the wave equation ([T]). They are determined 
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Thus the spatial dependence of all the optical parameters comes solely from the axial density 
profile of the condensate. In the following section we shall exploit the fact that it is only 
the density profile that determines the local optical properties of the condensate to develop 
an exact analytical solution of the wave Eq.([T]). 



III. ANALYTICAL AND SEMI- ANALYTICAL METHODS FOR ULTRASLOW 
PULSE PROPAGATION 

n 

For a uniform density medium the Eq{T] can be solved analytically [13]. The complex 
envelope of the incident wave is Gaussian pulse E{0,t) = e~*^/^o. The initial pulse, af- 
ter propagating in the medium of length L, is then found to be delayed with respect to 
a reference pulse propagating in vacuum by = L/vg. The final width of the pulse is 
r(L) = To-y/l + {L/ zqY with tq is the initial temporal width of the pulse, and zq = — tttq /62- 
For L ^ zq me get r(L) = \b2\L/TiTQ. Experimentally measured group velocity is de- 
fined by Vg = L/td, where the effective axial length of the medium is evaluated by 
L = [(47r/iV)/~ rdrC dzz' p{r, z)f' . 

For a non-uniform medium, the wave equation ([T]) with spatial dependent coefficients 
can be solved by using some basic differential equation solving methods. In the typical 
slow light experiments a small pin hole is introduced to couple incoming light with the 
condensate guide and axial propagation is enforced. In the subsequent discussions paraxial 
effects shall be ignored and strictly one dimensional propagation will be considered by taking 
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r = (x, z) = (0, 0, z). Fourier transforming Eq. ([T]) from t space to w space gives 

iti? 



^ + a{z)£ 
dz Vg{z} 



£ - w^ib2iz)£ = 



(7) 



where £ = £{z,w) = (l/v27r) E{z,t) exp {iwt)dt. Solution for the equation in w space 
becomes 



£{z, w) = £o exp 



(8) 



(^iw'^b2{z') + iw/vg{z') —a{z')^dz' 

where £o = (l/v^) JZo exp {-{t - tof/2T^) exp (iwt)dt = tq exp (-wVqVS). Transform- 
ing back to t space we find 



E{z,t) 



^0 



exp 



4/i(^) 



(9) 



where 



/2W 



!o _ 

t-to 



b2{z')dz' 
1 



(10) 
(11) 



Writing the trap potential with a variable curvature parameter k, such that V = nz"^, den- 
sity profile of ID Bose-Einstein condensate becomes p{z) = [(/i — V{z))/Uo]Q{fi — kz"^) + 
(73/2 (exp — /i)))/Ay. To calculate the pulse envelope analytically, it is necessary to 

be able to evaluate the single integral 



N{z) 



p{z)dz =N^{z)+Nt{z), 



(12) 



with Nq{z) = /—-pQ{z)dz and Nt{z) = fl^pT{z')dz'. Nq makes the dominant contri- 
bution in the condensate region. Let us define an auxiliary function for notational simplicity 
such that 
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where 0(x) = l/(y^) J^^e "^^^"^du is the normal cumulative distribution function which 
can be expressed in terms a tabulated special function, error function as 0(x) = 1/2(1 + 
erf(x/V2)). It is now possible to express the results of the integral as follows 

F{z) z < -\zo\ 



^ -F{z) z > \zo\ 



(14) 



Here \z()\ = J^/k. Finally we can rewrite optical pulse parameters in terms of N{z) so that 



This completes our analytic exact solution. Though it contains an infinite series, not all 
terms would be of significance at condensate temperatures of interest. Furthermore one can 
still make the result of more practical value by noting that it contains a special function. 
In addition to its series and asymptotic expansions, there are elementary, Gaussian-like 
functions that can be fit to the error function. Thus these facts encourage us to look for a 
semi-analytical method in which we fit polynomials or Gaussian to the N{z). Equivalently 
and more simply then one can make such fits to the density profile of the condensate. Its 
exact form is not essential as it only appears as an integrand. Due to the approximate fitting 
involved, this method would be a semi-analytical method to determine the pulse envelope. 
After quickly reviewing standard numerical solvers of wave equation such as Cranck-Nicolson 
and spectral methods, we shall test the semi-analytical method against them. 

IV. TYPICAL NUMERICAL METHODS FOR PULSE PROPAGATION 
A. Crank-Nicolson Method 

A dimensionless form of Eq. ([1]) can be solved via finite difference Crank-Nicolson space 
marching scheme. The Crank-Nicolson scheme is less stable but more accurate than the 
uUy implicit method; it takes the average between the implicit and the explicit schemes 
14| . Discritization is performed as follows, with i,j being the space and time grid variables 




(15) 
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If we plug them into the wave equation ([T]), we get 



f 1 ibo 



2vgAt 2 (At) 
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By adding the boundary conditions (in our work set as zero: Eq = E\^_i = for all "z") we 
obtain set of N linear equations with N unknowns, which have to be solved simultaneously 
for every space step i where the vector {Eq} is defined by initial conditions. Discrete 
equations in matrix form are solved using Thomas algorithm [l^ which is a fast Gaussian 
elimination method for tridiagonal matrices. 



B. Pseudo-spectral Method 

Instead of doing a finite difference approximation in time, we can expand the function 
E{z,t) in spectral series at a given position for all time values for better approximation of 
the time derivative. The initial value can be used to determine the coefficients of spectral 
series [isl. 



1^. An appropriate spectral series can be Fourier series. The reason we choose 



Fourier series instead of polynomials or any other series that the derivative of the Fourier 
series is just multiplication of Fourier series with a pre-defined vector and also fast Fourier 
algorithm makes it faster to compute the Fourier coefficients. 

The initial function is divided into N points (Ej) and discrete Fourier transform is applied 
in the interval -N/2 <k< N/2 - 1. 

1 N-l 

Eu = HEA = ^ E E.e-^'-'^/'' (20) 
i=o 

for all k. We get the following ordinary differential equation; 

^^ = -(a + |-a,F)B(^,^ (21) 

We solve this ordinary differential equation for discrete space intervals assuming that for 
each interval coefficients are constants; 

E{zo + Az, k) = E{z = zo, k) ■ ^^-(-"-f +^"2^^) (22) 
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V. RESULTS AND DISCUSSION 



In our numerical calculations, we specifically consider a gas of = 8.3 x 10^ ^^Na atoms 
for which M = 23 amu, Aq = 589 nm, 7 = 27r x 10.01 MHz, = O.57, = 27r x 10^ 
Hz, and = 2.75 nm. For the parameters of the trapping potential, we take Uj. = 271 x 69 
Hz and u;^ = 27r x 21 Hz as in Ref. The coupling field Rabi frequency is taken to be 
Qc = 0.567 [!]• Critical temperature for Bose-Einstein condensation of such a gas is found 
to be Tc = 424 nK. 

Fitting a polynomial of degree 22, pulse envelope is calculated using the analytical for- 
mulae. To illustrate the success of the fit we present the absorption coefficient in Fig.([T]) at 
temperature T = 42 nK. Solid line is the exact analytical absorption coefficient while the dot 
line is the semi-analytical result obtained after the polynomial fit. Similar behavior occurs 
for 62 and 1/vg. With the polynomial expressions of the optical parameters, the integrals 
required for pulse envelope functions such as fi^2, or N{z), are evaluated quickly and simply. 
Contour plots of the propagating pulse are shown in Figs. [SHHl We assume a Gaussian pulse 
with unit amplitude of the form exp (— (t — to)^/2Tg ) at initial time to, , where tq is the pulse 
width. When the optical pulse enters the condensate region, its group speed dramatically 
reduces under KIT conditions, as shown in Fig. [2l The optical pulse rapidly assumes its 
high speeds again following the passage to the thermal background cloud. We assumed 
the optical pulse is propagating in vacuum before and after the thermal component of the 
ultracold atomic system. In experiment, group velocity is measured in terms of time delay 
of the pulse with respect to a reference pulse which propagates in vacuum over the same 
distance with the atomic medium. Broadening of the pulse after leaving the condensate is 
visibly seen in the figure as it gets about almost two times broader. 

Similar behavior of the pulse but with a significant difference regarding the pulse width 
can be seen in the Fig. [31 The only change in the parameters used in Fig. [2] is that now 
Qc = 1.57. In that case the broadening and absorption becomes negligible while the pulse 
gets faster. These results, in particular the role of the intensity of the control field on the 
dispersion and loss management have already been discussed before in Ref. |9|. Here, the 
semi-analytical method is shown to reproduce them efficiently. 

To test the semi-analytical method with polynomial fitting against standard numerical 
solvers of the wave equation, a dimensionless form of Eq. is also solved via finite difference 
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Crank-Nicolson (C-N) space marching scheme and pseudo-spectral method (P-S). Final pulse 
width and final amplitude determined by these methods are listed in Tab.([T]). A microsecond 
pulse broadens by a factor of approximately ~ 1.7 according both numerical and semi- 
analytical methods as seen in Tab. ([I]). Similar agreement of the semi-analytical method 
with the numerical results are found for the absorption loss. 

As a further test on our numerical methods as well, we have compared the results of 
the Crank-Nicolson and pseudo-spectral codes against the exact analytical solution for a 
uniform density condensate cloud of p = 1.56 x 10^° 1/m^. The coefficients in the Eq. 
([T]) in that case are found to be Vg = 1.5 m/s, a = 2.1 x 10^ l/m, and 62 = 3.39 x 10~^ 
s^/m. Initial value for pulse width is 1 /isec. The final pulse width determined according 
to the Crank-Nicolson, pseudo-spectral, and completely analytic methods are 3.3863 /isec, 
3.3861 /isec, and 3.3869 /isec respectively. 

We can also use Gaussian fitting functions instead of polynomials in order to get more ex- 
plicit and compact expressions. Fitting a Gaussian to the density for the same experimental 
parameters, we have found the optical parameters become 



Reliability of the Gaussian fitting in the semi-analytical method is tested against spectral 
method and Crank-Nicolson method as summarized in Tab. flIT|) and Tab. fillip . Results ob- 
tained with numerical or semi-analytical Gaussian fit methods are found to be in good 
agreement among themselves. Furthermore, Gaussian fit method gives the similar results 
obtained with the polynomial fit. This should be the case as the density only enters as an 
integrand and the exact shape of the density should not be essential to determine the optical 
pulse properties. Gaussian fit method gives the final amplitude and width to be 0.475 and 
1.7524 yusec, respectively. In these tables we have seen that the agreement between the semi- 
analytical method and the numerical methods seems to improve as the grid made finer and 
finer. This may suggest that the semi-analytical method is more accurate than the standard 
numerical solvers. However we are unable to give more rigorous proof for that statement 
than these numerical tests. 




(23) 
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Having analytical pulse envelope expressions or its reliable and compact semi-analytical 
form allows us to vary controllable experimental parameters easily to optimize pulse width, 
group velocity and absorption loss of the pulse. To illustrate such optimum dispersion 
management, we choose to examine the effect of the trap curvature to illustrate. Effect 
of the control field intensity has been numerically investigated earlier ^]. Considering a 
quadratic trap as = nz^, trap curvature is translated to the curvature of the dielectric 
function through the density dependent EIT susceptibihty. Typical results obtained by the 
Gaussian fitting semi-analytical method are listed in listed in Tab. (]IV|) . The parameters are 
the same with those of FiglHand k is in units of kg/sec^. It is seen that pulse shape is better 
preserved at larger k. 

VI. CONCLUSION 

We have discussed analytical and semi-analytical solution of the one -dimensional wave 
equation which governs the propagation of an ultraslow optical pulse in a dispersive inho- 
mogeneous atomic condensate. Ignoring paraxial effects, slowly varying pulse envelope wave 
equation is solved by using Fourier transformation technique and by using special functions, 
in particular, normal cumulative distribution function which is related to the error function. 
It has been argued that, as the exact density profile only enters as an integrand to the ex- 
pressions, simpler polynomial and Gaussian fits can be made for more compact and simpler 
expressions. Such a semi-analytical method is heavily tested against standard numerical 
solvers of the wave equation, namely, Cranck-Nicolson and pseudo-spectral methods. The 
results confirm the reliability of the compact expressions obtained under semi-analytical 
method. As an illustration of the efficiency of analytical expressions, the role of trap curva- 
ture on the dispersion management has been investigated. A quick calculation shows that 
pulse shape is more preserved in traps with higher curvatures. Other experimentally con- 
trollable parameters can be similarly studied for optimized design of optical traps and EIT 
conditions to efficient storage of coherent optical information or to engineer ultraslow pulse 
shapes. 
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List of Figure Captions 

Fig. 1. Solid line shows the position dependence of the absorption coefficient a along 
the 2;-axis. The dot line is polynomial fitting for the absorption coefficient and the degree of 
the polynomial is 22. The ultracold atomic system of ^^Na = 8.3 x 10^ atoms at T = 42 
nK under EIT scheme. The parameters used are M = 23 amu, = 2.75 nm, Aq = 589 nm, 
7 = 27r X 10.01 MHz, = 0.57, = 0.567, Ta = 27r x 10^ Hz. 

Fig. 2. Contour graph of the propagation of a microsecond pulse through the interacting 
BEC by semi-analytical method. Time if) is scaled by 0.22 /is. and position i^z) is scaled by 
1 yum. The parameters are the same with those of Fig{T] 

Fig. 3. Contour graph of the propagation of a microsecond pulse through the interacting 
BEC by semi-analytical method. Time (t) is scaled by 0.22 /isec. and position i^z) is by 1 /im. 
The parameters are the same with those of Fig{T] except for fic = 1-57. 
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List of Table 



TABLE I: Comparison of Crank-Nicolson (C-N), pseudo-spectral (P-S), and semi-analytical (S-A) 
method. Propagation of optical pulse with initial pulse width 1/xsec and initial amplitude 1. 



C-N 


P-S 


S-A 


final pulse width 1.7311 //sec 


1.7305 iisec 


1.7309 iLtsec 


final amplitude 0.4755 


0.4754 


0.4758 



TABLE II: Ratio of the results for pseudo-spectral method to the semi-analytical values for the 
given position X time grid dimensions. 



grid dimensions 


final amplitude 


final width 


29x29 


0.998 


0.9993 


210 X 210 


0.9993 


0.9997 


212 X 210 


0.9995 


0.9998 


212 X 212 


0.9998 


0.99991 



TABLE III: Ratio of the results for Crank-Nicolson method to the semi-analytical values for the 
given positionxtime grid dimensions. 

grid dimensions final amplitude final width 

29 X 29 0.9998 0.9994 

2i2 X 2io 0.99988 0.99991 

2i2 X 2i2 0.99998 0.99994 
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TABLE IV: Dispersive propagation of optical pulse as a function of k. 



K 


final amplitude final width (/xsec) 


2 X 10-^2 


0.4113 


1.8916 


0.5 X 10-21 


0.5250 


1.6247 


1 X 10-21 


0.6121 


1.4612 


0.5 X 10-20 


0.7629 


1.2422 


1 X 10-20 


0.8212 


1.1741 
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FIG. 1: Solid line shows the position dependence of the absorption coefficient a along the z-axis. 
The dot line is polynomial fitting for the absorption coefficient and the degree of the polynomial 
is 22. The ultracold atomic system of ^^Na iV = 8.3 x 10^ atoms at T = 42 nK under BIT scheme. 
The parameters used are M = 23 amu, = 2.75 nm, Aq = 589 nm, 7 = 27r x 10.01 MHz, Ts = 0.5j, 
= 0.567, Ta = 27r X 10^ Hz. 
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FIG. 2: Contour graph of the propagation of a microsecond pulse through the interacting BEC 
by semi-analytical method. Time (t) is scaled by 0.22 ^us. and position (z) is scaled by 1 /xm. The 
parameters are the same with those of Fig(TJ 
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FIG. 3: Contour graph of the propagation of a microsecond pulse through the interacting BEC 
by semi-analytical method. Time (t) is scaled by 0.22 fisec. and position (z) is by 1 /im. The 
parameters are the same with those of Fig(T] except for Qc = I.57. 
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